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Abstract 

Multistationary chemical reaction networks are of interest to scientists and 
mathematicians alike. While some criteria for multistationarity exist, obtain¬ 
ing explicit reaction rates and steady states that exhibit multistationarity for a 
given network—in order to check nondegeneracy or determine stability of the 
steady states, for instance—is nontrivial. Nonetheless, we accomplish this task 
for a certain family of sequestration networks. Additionally, our results allow 
us to prove the existence of nondegenerate steady states for some of these se¬ 
questration networks, thereby resolving a subcase of a conjecture of Joshi and 
Shiu. Our work relies on the determinant optimization method, developed by 
Craciun and Feinberg, for asserting that certain networks are multistationary. 
More precisely, we implement the construction of reaction rates and multiple 
steady states which appears in the proofs that underlie their method. Further¬ 
more, we describe in detail the steps of this construction so that other researchers 
can more easily obtain, as we did, multistationary rates and steady states. 

Keywords: Mass-action kinetics, chemical reaction networks, 
multistationarity, determinant optimization method, steady states, degeneracy 


1. Introduction 

Although many dynamical systems arising in applications exhibit bistability, 
there is no complete characterization of such systems. Even for the subclass 
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of chemical kinetics systems and even under the assumption of mass-action 
kinetics, which is the focus of this work, the problem is difficult. 

Here we consider the simpler, yet still challenging, question: which chemical 
reaction networks are multistationary, i.e. which have the capacity to exhibit 
two or more steady-state concentrations with the same reaction rates? Mathe¬ 
matically, this asks: among certain parametrized families of polynomial systems, 
which admit multiple positive roots? Therefore, this is a real algebraic geometry 
problem, and we do not expect an easy answer in general. 

The first partial answers to this question are due to Feinberg, Horn, and 
Jackson in the 1970s. Their results in chemical reaction network theory mm 
(specifically, deficiency theory [3]) can preclude or guarantee multistationarity 
for certain classes of networks. For a survey of these and other methods, see [4]. 

Our work pertains to two related results: (1) a method for “lifting” multiple 
steady states from small networks to larger ones, and (2) the so-called determi¬ 
nant optimization method for certifying that a given network is multistationary. 

The lifting result, stated informally, is as follows: if a chemical reaction 
network contains an “embedded” network that is multistationary, then the entire 
reaction network also is multistationary under certain hypotheses [5] . Therefore 
we are interested in cataloguing the multistationary networks which contain no 
embedded multistationary networks, because all larger multistationary networks 
contain at least one embedded multistationary subnetwork from the catalogue. 

As a step toward such a catalogue, Joshi and Shiu identified a certain infinite 
family of chemical reaction networks Km,n to be of particular interest among all 
networks that include inflow and outflow reactions 0. This family is minimal, 
in that it has no embedded subnetworks (with inflow and outflow reactions) 
that exhibit multistationarity. To analyze these networks, Joshi and Shiu used 
the second method for analyzing multistationarity mentioned above. 

Developed by Craciun and Feinberg, the determinant optimization method 
can assert that a network is multistationary 012]; as such, it is a partial converse 
to their results on “injective” reaction networks which guarantee that a network 
is not multistationary. This topic of injectivity has seen much interest in recent 
years (see 0IHIIS] and the references therein); however, the determinant opti¬ 
mization method has garnered comparatively little attention. The only related 
results that we are aware of are due to Banaji and Pantea [8], Feliu [10], and 
Miiller et al. [Sj. 

Using the determinant optimization method, Joshi and Shiu proved that 
Km.n is multistationary for all integers m > 2 and odd integers n > 3 (Ij. 
Furthermore, they conjectured that these networks can exhibit multiple nonde- 
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generate steady states. (It is not guaranteed that the determinant optimization 
method produces nondegenerate steady states; we show this for the first time 
in Remark 3.11| ) The significance of the conjecture is that if it is true, then 
K 2 ,n would be the first example of an infinite family of chemical reaction net¬ 
works with inflow and outflow reactions and at-most-bimolecular reactants and 
products—that is, minimal with respect to the embedding relation among all 
such networks which have the capacity to exhibit multiple nondegenerate steady 
states. Nondegeneracy is important because results that “lift” multiple steady 
states from embedded subnetworks or other typically smaller networks require 
the steady states to be nondegenerate; a summary of such results appears in [H 
§4]. Also, because trimolecular reactants/products are rather uncommon in 
chemistry and K 2 ^n is at most himolecular, this family of networks is of partic¬ 
ular interest in chemical applications. 

In the current work, we resolve the conjecture for the case n = 3 and all 
m > 2; in other words, we prove that Km ,3 has the capacity to admit multiple 
nondegenerate steady states for all m > 2 (Theorem |4.5[ ). To accomplish this, 
we need information beyond the mere existence of multiple steady states; we also 
need precise values (or at least estimates) for the rates and steady states. By 
applying the proofs underlying the determinant optimization method in Craciun 
and Feinberg’s work to the networks Km,n, we obtain (via standard methods for 
analyzing recurrence relations) explicit closed forms for multistationary rates 
and steady states. Then we use these closed forms to verify that the steady 
states are nondegenerate for small values of m and n. 

Finally, recognizing the usefulness of generating closed forms (or at least 
estimatetQ) for rates and steady states for any reaction network that satisfies 
the hypotheses of the determinant optimization method, in Section we outline 
the steps of the method with enough generality to be used in other contexts. 
These steps are present in Craciun and Feinberg’s work but are spread out 
over several proofs, so our contribution here is to reorganize the method into a 
concise procedure. 

An outline of our work is as follows. Section introduces chemical systems 
and the main conjecture. Section describes the determinant optimization 
method in detail. We use this method to resolve some cases of the main conjec¬ 
ture in Sections 1^ andFinally, a discussion appears in Section]^ 

Notation. We denote the positive real numbers by M_|_ := {z G M | x > 0}, 


^For general networks, the determinant optimization method need not yield closed forms 
for the rates and steady states, but one can nonetheless obtain estimates. 
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and the standard inner product in M" by The entry of a vector x is 

denoted Xi. 


2 . Background 

This section introduces chemical reaction networks, their corresponding mass- 
action kinetics systems, and the main object of our paper: the sequestration 
network Km,n- 

2.1. Mass-action kinetics systems 

Definition 2 . 1 . A chemical reaction network G = {S,C,TZ] consists of three 
finite sets: 

1. a set o/species S = {Xi,X 2 ,... ,-^s}; 

2. a set C of complexes, which are non-negative integer linear combinations 
of the species, and 

3. a set TZ C C X C o/reactions. 

Example 2.2. The following chemical reaction network: 

Xi+X2-^ X3 
X2 — y X\ X4 , 


is entirely defined by: 

1. the set of species S = {Xi, X 2 , X 3 , X 4 }, 

2. the set of complexes C = {Xi -|- X 2 , X 3 , X 2 , Xi -\- X 4 }, and 

3. the set of reactions TZ = {(^1 -I- X 2 , X 3 ), (X 2 , Xi -\- X 4 )}. 

Any reaction network G = {S,C,TZ} is contained in the fully open extension 
network G obtained by including all inflow and outflow reactions: 

G := {5, C U 5 U {0} , 7^U {X, o . (1) 

In other words, the fully open extension of any network is obtained by adding 
the reactions X —> 0 (inflow) and 0 —> X (outflow) for all X G S. 

As all the reactions take place, the concentrations of each of the species 
will change. We make use of mass-action kinetics to define a system of ordi¬ 
nary differential equations that describes, for each species, how its concentration 
changes as a function of time. This ODE system is described by the stoichio¬ 
metric matrix T and the reactant vector R{x), which is a vector-valued function 
of the vector of species concentrations x. 
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Definition 2.3. Let G = {S,C,TZ} be a network, and let {yi —>■ —>■ 

j/ 2 ) • • •, y\'R\ v'\n\^ ordering of the reactions. 

1. The reaction vector of the reaction yi —>■ y^ is the vector y'i — yi, viewed in 

Note that yi y[ is a slight abuse of notation, used to denote the 
reaction yi • X y[ ■ X where X is the vector of all species. Explicitly, 
the vectors yi and y'i only contain species coefficients. 

2. The stoichiometric matrix of G is the |5| x \R.\ matrix T whose column 
is the reaction vector of yk ^ u'k- 

3. The reactant vector R{x) is the vector of length \R.\ whose entry is the 
(monomial) product: 


rkX'i 


(Vk)lXVk)2 


(yk)i. 

151 


where rk G R+ is the reaction rate of the kf^ reaction. 


Definition 2.4. The mass-action kinetics system of a network G = {S,C,TZ] 

17? I 

and a vector of reaction rates (jk) G is defined by the following system of 

ordinary differential equations: 

g = r.i^(x). (2) 

Example 2.5. For the following network: 


A + 2B^2A , 


r = 


1 

-2 


and R{x) = (rWAx"^), so the mass-action kinetics system (2) is 


dXA 

dt 

dxB 

dt 


r • R{x) 


rxAX% 

—2rxAX% 


An important characteristic of mass-action kinetics systems is that they may 
or may not have the capacity to admit (positive) steady states: 

ISI 

Definition 2.6. A positive steady state is a vector X* G such that r • 
i?(x*) = 0. A steady state x* is nondegenerate */ Im(d/(x*)|ii„(r)) = Im(r), 
where dffx.*) denotes the Jacobian matrix of the mass-action kinetics system at 

X*. 
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Definition 2.7. A network that includes all inflow and outflow reactions is mul- 
tistationarj0 if there exist two distinct concentration vectors x*, and positive 
reaction rates such that F • i?(x*) = F • R{x^) = 0. 

2.2. The sequestration network Km,n 

The main object of our paper is the fully open extension of the network 

Definition 2.8. For positive integers n > 2 and m > 2, the sequestration 
network Km,n is: 

Xi+X2^0 (3) 

X2+X3^0 


Xn-l + Xn —>■ 0 

Xi —^ uiXji . 


Km.n is the fully open extension of Km,n, obtained by adjoining all inflow and 
outflow reactions, as in Q. 

Schlosser and Feinberg analyzed variations of iF 2 ,n El Table 1], as did 
Craciun and Feinberg 0 Table 1.1]. Joshi and Shiu introduced the version of 
the sequestration networks in Definition |2.8[ and proved that some of them are 
multistationary: 


Proposition 2.9. Lemma 6.9] For positive integers m>2 and n > 3, if n 
is odd, then Km,n admits multiple positive steady states. 


For m = 1 or n even, the network ^ is “injective” and therefore not 
multistationary §6]. Joshi and Shiu conjectured that Proposition 2.9 extends 
as follows: 


Conjecture 2.10. For positive integers m > 2 and n > 3, if n is odd, then 
Km,n admits multiple nondegenerate steady states. 


To resolve Conjecture 2.10 we must show that Im{df{x*)) = Im(T) for two 
distinct positive steady states x*. We will see in Q below that F is full rank, 
so we need only show that det{df{x*)) 0 for two positive steady states x*. 


^The focus of this work is on certain networks Kmn, that include all flow reactions. For net- 
works that do not include all flow reactions, the definition of multistationary must incorporate 
the conservation relations in the network, if any. 
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Remark 2.11. As mentioned in the introduction, networks in which all reac¬ 
tants and products are at most bimolecular— that is, each complex has the form 
0, X, X Y, or 2X — are the norm in chemistry. This is the case for the 
networks 1 ^ 2 ,n; so that the internal reaction is Xi —> 2X„. 


We end this section by displaying the matrices that define the mass-action 
kinetics system ([^ defined by Km,n- We order the reactions as follows: first, we 
enumerate the n internal (or true) reactions listed in (|^ (so, the first reaction 
is Xi X 2 —t 0 , and so on), next are the n outflow reactions (so, the (n -|- l)st 
reaction is Xi —> 0 , and so on), and then we have the n inflow reactions (so, the 
(2n -I- l)st reaction is 0 —>■ Xi, and so on). We will refer to the sets of internal 
(true), outflow, and inflow reactions as TZt, 'R-o, and TZi, respectively. 

The stoichiometric matrix for Km,n is: 


-1 

0 

0 .. 

. 0 

-1 



-1 

-1 

0 .. 

. 0 

0 



0 

-1 

-1 ■■ 



-in 

jn 

0 

0 

-1 ■■ 








. -1 

0 



0 

0 

0 .. 

. -1 

m 




(4) 


where is the n x n identity matrix. The reactant vector is: 


riXiX2 

r2X2X3 


R{k) 


Xn—lXn—lXn 

XnXl 

Tn+lXl 

rn+2X2 


X2nXn 

r2n-l-l 


ran 


where the ri G K+ are the reaction rates and each Xi G M+ is the concentration 
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of each species Xi. The mass-action ODEs (2.41 are: 


±1 =- riXiX2 - TnXl - Xn+lXi + r2n+l 

Xi =- Ti^iXi-iXi - TiXiXi+i - Xn+iXi -I- r2n+i for 2 < i < n - 1 

in =- Xn-lXn-lXn + TUT^Xi - r2nXn + T’3« • 


Thus the Jacobian matrix, df{x), is the following (n x n)-matrix 


'—riX2 -Vn- r„+i 

-riXi 

0 


0 

0 

-r-iX2 

-TiXi - r2X3 - r„ + 2 

-r2X2 




0 

— r2X3 

-1-2X2 - T3X4. - Tn+Z 


0 

0 


0 

-r3Xi 


-rn-2Xn-2 

0 

0 

mrn 

0 

0 


-rn-2Xrt-2 “ Tn-lX^ - r2n-l 

— Tn-lXn 

—Tn-lXn-l 

— rn-lXn — 1 — T2n. 


(5) 


3. Constructing multiple steady states via the determinant optimiza¬ 
tion method 

The determinant optimization metho^was developed by Craciun and Feinberg to 
show that certain chemical reaction networks are mnltistationary [6]. More precisely, 
the method gnarantees that some networks (such as those that satisfy the ‘Input’ 
conditions below) are necessarily multistationary. For instance, Joshi and Shiu showed 
that the networks Km,n (for m > 2 and odd n > 3) satisfy the ‘Input’ conditions, and 
thus concluded these networks are multistationary (Proposition |2.9| . 

In fact, the determinant optimization method also applies to some networks that 
do not satisfy the ‘Input’ conditions. To determine if this is the case for a given 
network, one must check whether a certain optimization problem has a solution (see 


Remark 3.41. If so, then the method guarantees that the network is multistationary. 


However, in many applications, it is useful not only to know that a network is 
multistationary but also to have explicit steady-state concentrations and reaction rates 
that are witnesses to multistationarity. For instance, here we would like to determine 
whether the steady states are degenerate, whereas in other settings one might like to 
perform stability analysis. 

Fortunately, the proofs in [6l §4] that underlie the determinant optimization method 
are constructive, up to one use of the Intermediate Value Theorem, so one can generate 
or at least approximate steady states and rates. This section describes the step-by- 
step procedure to do this; following our steps is easier than (although equivalent to) 
“backtracking” through the proofs. That is, our contribution here is to re-package 
the determinant optimization method into a constructive algorithm. We will see that 


^Related techniques for establishing multistationarity appear in work of Banaji and Pan¬ 
tea [S] §4], Feliu |10l §2], and Muller et al. [O] §3.2]. 









for some networks, such as Km,n, the method constructs closed forms for the steady 
states and rates. 


Determinant optimization method (constructive version) 

Input: Any chemical reaction network G with n — |<S| species that contains all n 
inflow reactions such that there exist n reactions yi —>■ t/i, 1/2 —>■ 2/L • ■ • > Vn y'n 
among the internal (true) and outflow reactions TZt U TZo of G for which 

(I) det(?/i,i/ 2 ,...,j/n) •det((j/i -y'i),{y 2 - t/ 2 ), (i/n - y'n)) < 0, and 
(II) there exists a vector ij G R" such that E(Li?)i(j/i — y'i) G R+ = r1^'. 

Output: A certificate of multistationarity of G: (approximations of) a positive re¬ 
action rate vector {ry_,yi) £ r!^' and two positive concentration vectors x* and 
which are both steady states of the mass-action system defined by G and {Vy^yi). 
Steps: Described below. 


With an eye toward resolving Conjecture 


2.10 


K„ 


will be our ongoing example. 


Example 3.1. For Km,n (with m >2 and n > 3 odd), hypothesis (II) is satisfied by 
the vector y = (1,1,... ,l,m + 1,1) ^ Lemma 6.9]. Hypothesis (I) was proven in 
Lemma 6.7], where the n reactions are precisely the n internal (true) reactions ©• 
Conveniently, these are the reactions labeled yi —^ y) of the sequestration network, for 
1 < i < n, so our notation for the first n reactions—as well as the use of n for the 
number of species—matches that of the determinant optimization method. 


The steps below involve a certain linear transformation specifically, for rj £ 
j^TZj-uTZo^ the linear transformation Ty : R''^' —^ R''®' is defined by: 

ny^y'(y-5){y-y')- ( 6 ) 

y^y' ClTItOTZo 

Equivalently, the matrix representation of Ty is d(—/)(1, !,...,!) where the rates are 
given by ri = ry. In other words, this matrix is the Jacobian matrix of the mass- 
action system (2.41 defined by the internal and outflow reaction rates rj (and any 
choice of inflows: they do not appear in the Jacobian matrix) at the concentration 
vector (1,1,..., 1). 


Example 3.2. For Km,n, the matrix representation of Ty is: 

m +ri„+ rj„+i Tji 0 

rji yi +rj 2 + rjn +2 r ]2 

0 rj 2 ri2 + rj3 + rjn+3 


V3 


0 


0 

-mrjri 


rjTi-2 -(- rjn-l -(- rj2n-l Vn-l 

rjn-1 rjn-l + rj2n .. 


(7) 

From the Jacobian matrix ([^, it is clear that this matrix ([^ equals d{—f){l, !,...,!), 
where the reaction rates are given by r; = ry. 
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The first step is to construct a (strictly positive) vector rj € indexed by 

all internal (true) and outflow reactions, such that 
(I) det(r,j-) < 0, and 

(11) E 

y^y' eTlT^T^O 

Craciun and Feinberg proved that these conditions (I) and (II) are satisfied by a vector 
ri~ of the following form: 


% 


Xfjy^y' if y^y' ^{yi^y[\i^ H) 

I e else , 

where A is sufficiently large and e is sufficiently small [6l proof of Theorem 4.2]. 
Example 3.3. For Km,,n (with m > 2 and n > 3 odd), we define rj~ as follows: 

A ifl<i<n — 2ori = n 

Vi = \ [m + 1)A if i = n — I 

^e if n + 1 < i < 2n . 


Remark 3.4 (Stronger versions of the determinant optimization method). This sec¬ 
tion describes the simplest version of the determinant optimization method. In fact, 
even if a network does not satisfy the hypotheses in the input given above, the method 
can still apply: Remark f.l] describes, in this setting, how to implement the above 

first step, i.e. how to test whether a suitable r]~ exists, as an optimization problem. 
Specifically, this is a polynomial optimization problem with linear constraints over a 
compact set. Therefore, one can use any applicable optimization method. Additionally, 
see Remark \3.5\ for how one can begin the algorithm at the second step. 

The second step is to construct a (strictly positive) vector rp G which: 

(I') det(r^o) = 0, and 

(11') E vl^yPy-y')(^^'+'- 

y^y' GIZtUTZo 

Craciun and Feinberg proved that this can be accomplished as follows [51 proof of 
Theorem 4.1]. First, construct an 77 + € (Jet(r^+) > 0; do this by 

assigning a large value to outflow reactions and a small value to internal reactions: 


V 


+ 

y-i-y' 


A"*" if 1/ —^ 1/' G TZo 
A y ^ y' G TZt , 


where A^ > 0 is large and > 0 is small. 

Then, by interpolating between this vector and the vector rj~ from the pre¬ 
vious step, the Intermediate Value Theorem (plus the fact that the set of vectors 
satisfying condition (II) is convex) guarantees the existence of an if with the required 
properties. Moreover, such a suitable rf can be numerically approximated, and, with 
careful tracking of error, one can use this approximation to generate steady states and 
concentrations in the following steps. 
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Remark 3.5 (Interpretation of the second step and subsequent steps). What the 
second step does is to find reaction rates (given by rjo for the internal and outflow 
rates, and the vector in {II') for the inflow rates) at which the concentration vector 
(1,1 ,...,!) is a degenerate steady state: degeneracy is by {!'), and being a steady state 
comes from {II')- 

Equivalently, if {vy^yi) is any positive vector of reaction rates at which some con¬ 
centration vector c is a degenerate positive steady state, then the vector p G 


p|7?.tU7?.o 


defined coordinate-wise by rjy- 


'(? satisfies the second step. So, a reader 


who has already found a degenerate positive steady state of their system could start the 
determinant optimization method at the second step. In other words, one could begin 
applying the method by immediately searching for a suitable rp (without first generating 
r}~ and ). One strategy for doing this is described in Remark 
for Km,n beginning in Example 


3.1 


3.6 


which we employ 


In the next steps, the determinant optimization method constructs a certain vector 
S so that |5| is a suitable bifurcation parameter: for |5| small but positive, the degenerate 
steady state breaks into two nondegenerate steady states. 

Remark 3.6. Here is one strategy for constructing a suitable rf (without using rj~ 


and rj'^). First, identify (if possible) an p £ and a reaction yi 
internal (true) and outflow reactions such that: 


y'i among the 


Vv^y'{y-y') e 


|S| 


(a) E 

y^y' e(TiT'~nio)\{yi^y'i} 

(b) Vi — y'i £ R^fg (this holds, for instance, if y.. 


and 


y'i is an outflow reaction). 


One could see whether or rj~ might work (we use ri~ in Example 3.1 below). Then, 
define rj as follows: let the entry rP (corresponding to the same reaction) be free, and 
fix rp = 77“ for j 7^ i. Then, solve the (univariate polynomial) equation det(r^o) = 0. 
If there is a positive solution (for rP), then the resulting vector rp is positive, and {!') 
holds by construction. Furthermore, {II') holds because the sum in {II') is precisely the 
sum of a positive vector (namely, the sum in (a)) and a non-negative vector (namely, 
rp{yi — y'i)). However, rp is not guaranteed to be positive, so this strategy may fail. 

Example 3.7. For Km,n (with m > 2 and n = 3, 5, 7, 9,11 ), the following choice of 
rp satisfies the requirements of the second step: 

'a 

(m + 1)A 


Vi = 


( 8 ) 


ifl<i<n — 2 or i = n 
if i = n — 1 
ifn-\-l<i<2n — 1 
if i = 2n , 

where A > 0 and e > 0 are such that rpn is positiv^ (thus, all coordinates of rp are 


(m + l)(mA"+A^{m+l)"l„_2) \^ 

(a(™+2)+.)-i„_2-a2-i„_3 - + 1) 


■^We checked that such A,e exist for n = 3,5,7,9,11 (and m > 2), and we furthermore 
conjecture that for larger n, choosing A sufficiently large and e sufficiently small will suffice. 
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positive), and ~[i is i*^ principal minor of the matrix representation of T^o displayed 
below in ([^, i.e. ~li is the determinant of the ixi (tridiagonal) upper-left submatrix 
of (also, “lo := ^)- 

To show that this choice of 77 ° satisfies the two conditions of the second step, we 
first note that {II') is straightforward to verify. 

Satisfying {!') only requires rf to satisfy one (determinantal) equation, so allowing 
one free variable is sufficient. We choose ri 2 „ as this free variable, and we will recover 
the formula in (|^. Namely, we let all other coordinates rf) have the form given above 
(for 1 < i < 2n — 1), and then, recalling 0 ’ the matrix representation of T^o is: 


'2\ + e A 0 

A 2A + e A 0 

0 A 

: 0 


0 

0 


0 


0 : 
—mA 0 


0 


A “t" \{m “t" 1) “t" 6 
A(m + 1) 


A(m + 1) 
A(m + 1) + 


( 9 ) 


Expanding ^ along the bottom row, we obtain the determinant ofT^o: 
det(r^o) = (-l)"m(m + 1)A” - A^(m + + (A(m + 1) + , (10) 


where we recall that ~[i is i*^ principal minor of T^o. 

The determinant of tridiagonal matrices can be solved recursively flSI) . For our 
matrix, it is easy to verify the following recursion for i < n — 2, which is independent 
of m: 


~li+2 — (2A + — A ~li , 


with initial values ~lo := 1 and “li = 2A+e. Notice that ~l„_i must be treated separately, 
because the (n — 1 )®* row contains 77°_i, which is a function of m. Using standard 
methods we can get the generating function of the recurrence. 

“li = 2^IT(r * (“C2(C2 - Cl)* + Cl(C2 - Cl)* + C2(ci + C2)* + Ci (ci + C2)*^ , 


where ci = {e)^{e + 4A)5 and C 2 = e + 2A. Notice here that “li is always positive 
for sufficiently small e. A formula for “l„_i is given using the formula for tridiagonal 
matrices 

“In-l = “ln-2(A(m + 2) + e) — A^“1„_3 . 


With this recurrence solved, we set equation (_?0| ) to zero and then derive the explicit 
function for rffn from ,10) in terms of m, A, and e; this is the formula in ®. 


The third step is to construct a nonzero vector 5 € in the nullspace of T^o, i.e. 
such that T^o -(5 = 0. (Such a vector 5 exists because det(T^o) = 0.) 
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Example 3.8. For our example Km,n (with m > 2 and n > 3 odd), we claim that the 
vector 5 whose coordinates are defined as follows is in the nullspace of 


5i 

Sk = { -4_2 

-(A(m + 2) + e) 


rSn-2 


if k = l 

if 2<k<n — 1 
if k = n , 


( 11 ) 


A(m+1) 

where we introduce (5o = 0 for convenience in solving the recurrence relation, and 
7 ^ 0 is our free variable. Note that the last coordinate, 5n, has a different formula, 
because the {n — 1)®* row in used to define 5n contains terms dependent on m that 
do not satisfy the recurrence. 

To see that 5 is a nonzero vector in the nullspace o/T^o, notice that the conditions 
on 5 that state that its inner product with each of the the first n—l rows o/T^o coincide 
precisely with the n—l recurrences in the definition of S We claim that the last 

row ofTrj is linearly dependent on the other rows, and from this we will conclude that 
the last row of '^vo automatically has zero inner product with S. To see this, we recall 
that det T^o = 0 by construction, so we need only show that the first n—l rows of T^o 
are linearly independent. Assume, to the contrary, that there is a non-trivial linear 
combination of the first n — I'’* rows of T^o that adds to 0. Note that in the first 
n—l rows only the last one contains an entry in the last column, namely ri),_i. Since 


entire row. In the same manner, the corresponding scalar for the n — row would be 
equal to 0. Repeating the process shows that the only linear combination that adds to 
0 is the trivial one, thus, arriving at a contradiction. So, 5 is in the nullspace ofT^o. 

Again, by using standard techniques for analyzing recurrences, we find the gener¬ 
ating function for each of the first n—l entries of 5: 


Sk = (5iA 

for 1 < k < n — 1. 


(x/4Ae + £2 - (2A + e))*^ - (-y/dAe + 6^ - (2A + e))* 
2''A*V4Ae + e2 


The final step is to use the vectors and S (or, as we will see, a sufficiently scaled 
version of 5) from the previous two steps to construct a certificate of multistation- 
arity proof of Lemma 4.1]; namely, the internal (true) and outflow reaction rates 
are: 

(y’ 


Vy^y' 


for all y ^ y' G TZt U TZo 


'y^y e<yA)-i' 
the inflow reaction rates are the coordinates of the following vector: 

p|S| 

'ly^y' \y y J ^ 1 

y-fy' 

and the two steady states are: 


{ro^xj = yv^y'^y~y'^ ^ 


( 12 ) 


= (1,1,-,1) 


and 


= (e^^e^^...,e^|5l) . 
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Craciun and Feinberg showed that for sufficiently small scaling of 5, all inflow rates (121 
are positive [6]. 


Example 3.9. For Km,n, with m >2 and n = 3, 5, 7, 9,11, we checked that 5i = 1 
sujfices. Details for the n = 3 case are provided in Remark \4-^ 

Summarizing what we accomplished above, we have closed-form expressions for 
reaction rate constants and steady states that show that the sequestration network is 
multistationary: 


Theorem 3.10. Consider positive integers m > 2 and n G {3,5,9,11}. Let 5 G 
be as in 


11) with (5i = 1. Also, let rf be as is fslr] Then, for the following internal 


(true) and outflow reaction rates: 


and the following inflow reaction rates: 


Pi for all i G {1, 2,..., 2n} , 


r 2 n+i =ri+rn + r„+i 

r 2 n+i = ri-i +ri-\- r„+i for all 2 < i < n — 1 
ran = Vn-I + r 2 n - mr„ , 

the concentrations: 

X* = (1,1,...,1) and x'^ = (e‘*\ e'*") (13) 

both are positive steady states of the mass-action kinetics system defined by Km,n and 
the reaction rates rt above. 


Remark 3.11. One may wonder whether or not the determinant optimization method 
has the potential to create degenerate steady states. Indeed, if we could prove that 
this method always constructs nondegenerate steady states, then this would resolve 
Conjecture \2.10\ However, this is not the case. 

We determined this by analyzing K 2 ,a as in Theorem 


3.10 


By letting e be a free 

variable we compute the parametrized determinants det(d/(x*)) and det(d/(x’^)) as 
functions of e. Both functions are easily checked to be continuous for positive values 
of e, and from the graph (Figure^, we can see easily that there exist choices of e 
for which one of the two steady states is degenerate. More precisely det(d/(x*)) = 0 
for some e G (0.12,0.125) and det(ci/(x’^)) = 0 for some e G (0.240,0.241) and some 
eG (1.159,1.160). 


®In fact, this theorem will hold for any larger n for which the last coordinate of 77 ® as 
defined in can be made to be positive. 
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Figure 1: Graphs of det(c(f(x*)) solid and det(c(/'(x#)) dashed as functions of e 


4. Resolving Conjecture |2.10| for the n — 3 case 


Recall that Conjecture 2.10 asserts that Km,n admits multiple nondegenerate pos¬ 


itive steady states, for integers m > 2 and n > 3 with n odd. The main result of this 
section (Theorem |4.5| l resolves the conjecture when n = 3. To accomplish this, we first 
write down rate constants for this n = 3 case for which there are two steady states 
X* and x’^; these values were obtained by the determinant optimization method in 
the previous section for the general Km,n case (Proposition 4.11. We then resolve the 
conjecture for n = 3 by proving that x* and x’^ are nondegenerate. 


4.I. Reaction rate constants for which Km ,3 is multistationary 

Proposition |4.1| below specializes Theorem |3.10| to the n = 3 case. Following the 
description in SectionA = 1 and e = 0.1 will suffice, and then we obtain 


Vo 


^A, A(m + 1), A, e, £, 


-0.31m- 1.31 
2.1m-h 3.41 


T 


from the second step of the determinant optimization method. Next, in the third step, 
we find that the following vector spans the nullspace of T^o: 


5 = 


2.1m-f 3.41 
m -|- 1 



Thus, Theorem |3.10| specializes to: 
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Proposition 4.1. Consider any integer m > 2, and the following internal, outflow, 
and internal reaction rates: 


^4 = ^ « .06 


r2 


1.31 
1.31 - 

e + 1 _l 


rs = 


-.21 

e-2.1-1 


.24 


rs = 


.58 


re 


m —1.31 

2.1m+3.41 
e m+l 


r? = ri + rs + n ~ 2.29 rg, = n + r2 + rs 


rg = r2 + re — mrs 


(14) 


r/ien /or t/ie mass-action kinetics system defined by the fully open seguestration net- 

~ Jt f —2 1 2.1m+3.. 

work Km ,3 and the above rate constants ri, bothx* — (1,1,1) andx^ = (e, e ' ,e ”*+4 
are positive steady states. 


Note that in Proposition 4.1 only x.f, rg, re, rg, and rg depend on m. 


Remark 4.2. The only reaction rate in (141 that is not obviously positive is the inflow 
rate rg, so we verify it here: 


rg = r 2 + re — mrg > r 2 + 0 — m 


1 

e — 1 


> m — m 


1 

e — 1 


> 0 


where the second-to-last inequality follows from Lemma\4.3\ below. 


4-2. Bounding rates and steady states of Km ,3 

Here we give upper and lower bounds which we will use to prove that x* and 
are nondegenerate. The following bounds are on the third coordinate of x'^: 

e y+i > Xg = e "*+4 > e^' for all m > y > 0 . (15) 

_ 2.1m+3.41 

The first inequality in (151 follows from the easy fact that e "*+4 is a decreasing 
function when m > 0, and the second inequality is straightforward. 

The proofs of the following two upper/lower bounds are in Appendix A: 

Lemma 4.3 (Bounds on rg)- When A = 1 and t = 0.1, the rate constant rg defined 
in (14l satisfies the following inequalities for all m > 2; 

m+l > rg > m . 

Lemma 4.4 (Bounds on rg). When A = 1 and e = 0.1, the rate constant rg defined 
in (141 satisfies the following inequalities: 

0.14m > rg > 0.13m —0.5 , 

where the upper bound holds for m > 2, and the lower bound holds for m > 20. 


16 


















4-3. Proving nondegeneracy of steady states for the network Km ,3 
The main result of this section is: 


Theorem 4.5 (Resolution of Conjecture 2.10 when n = 3). For integers m > 2, the 
network Km ,3 has the capaeity to admit multiple nondegenerate positive steady states. 


We will prove Theorem 4.5 by showing that x* and in Proposition 


4.1 


nondegenerate, i.e. we must prove that the image of the 3x3 Jacobian matrix d/(x) 
at each of the steady states is equal to the image of the 3x9 matrix P. As stated 
earlier (after Conjecture 2.101, since P is full rank, our problem reduces to showing 
that det(ci/(x)) 7^ 0 for both steady states and for all integers m > 2. 

We begin by displaying the Jacobian matrix of Km,3- 


rf/(x) 


—r-iX2 — rs — Vi —ri*! 0 

—riX2 —rixi — r2X3 — rs —r2X2 

mrs —r2X3 —^ 2*2 — re 


Thus, our goal is to show that the following determinants (obtained by strategically 
cancelling and rearranging terms) are nonzero for all integers m > 2: 


Di := det(d/(x*)) = r 2 r\r 3 m — {r 2 + r(,){rir 3 + nn + r^rs-\-rsvs-\-nrs) 

— r2r6(ri + rs + r4) (16) 

D 2 ~ det{df(yi*)) = r2xf{{rixf + rs + r4,){r2xf) + rixfmrs) 

— {r 2 xf +re){rixf +r3 + r4)(rixf + r 2 xf + rs) 

+ {r 2 xf +re){rixfrixf) (17) 


From its graph|^ Di appears to be increasing quadratically as a function of m. So, 
to prove that Di is nonzero for integer values of m > 2, we will bound it from below by 
a quadratic function. Similarly, D 2 appears to be decreasing quadratically, so we will 
bound it from above by another quadratic function. Using these bounds, we will then 
conclude that Di and D 2 are strictly positive and negative (respectively) after certain 
cutoff points of m, effectively showing nondegeneracy of both steady states beyond the 
cutoffs. Finally, we will complete the proof by evaluating Di and D 2 at the remaining 
integers m between the 2 and the cutoff points to show that these values are nonzero. 


Proof of rfeeorem |J.5[ We generate our rates and concentrations as in Proposition [4^ 
Following the description immediately after Theorem |4.5| we need only show that 


Di 7^ 0 and D 2 7^ 0 for all integers m > 2. First, we bound Di by using its formula (161 
together with the bounds in Lemmas |4.3| and |4.4| 


Di > m^rira — ((m + 1) + 0.14m)(rir3 + rir4 + rirs + rsrs + r4r5) 
— (m + l)(0.14m)(ri + r3 + r4) . 


® Analogous graphs for larger n appear in Appendix B. 
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Next, estimating the remaining rates Vi, which are constants (recall equations (14l), 
by appropriate upper or lower bounds, we obtain: 


Di > m^(0.95) - ((m +1)+0.14m)1.61 - (m+l)(0.14m)2.29 
= 0.6294m^ - 2.156m-1.61 . 


It is easy to show that the quadratic function which bounds Di above is always positive 
for integers m > 4. So, Di > 0 for m > 4. Thus, it remains only to show that D\ ^ Q 
at m = 2, 3,4; indeed, those values are nonzero and are listed in Table 


m 

2 

3 

4 

5 

6 

7 

8 

9 

10 

Di = 

det(ci/{x’)) 

0.336 

2.784 

6.525 







D 2 = 

det(d/(x#)) 

-1.063 

-3.811 

-7.85 

-13.19 

-19.8 

-27.71 

-36.89 

-47.36 

-59.11 

m 

11 

12 

13 

14 

15 

16 

17 

18 

19 

D 2 = 

det(d/(x#)) 

-72.14 

-86.4 

-102 

-118.9 

-137.1 

-156.5 

-177.2 

-199.2 

-222.5 


Table 1: Determinants of the Jacobian matrices (16-17) at the two steady 
states X* and for the values of to > 2 before the proven bounds are valid. 
All of these determinants are nonzero, so the corresponding steady states are 
nondegenerate. 


Now we proceed to bound D2. Again, we use its formula 0 together with the 


bounds in (15) and Lemmas 4.3 and 


4.4 


D 2 < {m-\-l)x*{(rixf-\-r 3 + ri){(m-\-l)x*)+ rixfmri) 

— {mx 2 + (.13m — .5))(ria:f + rs + r4)(ria;f + mxf + rs) 

+ {{m+Vjx* + .14m)(ria;friif’) . 

Note that we used the lower bound on rg, so the above inequality holds for m > 20 
(and thus we will need to check the values of m between 2 and 19 separately). In the 
same manner as before, we approximate all of the constants appropriately for m > 20, 
and then simplify: 

D 2 < (m + l).13((.85)((m+ 1)8.7) +2.61m)-(.25m-.5)(.84)(4.72 + 8.16m) 

+ (.13(m+ 1) + .14m)(.91) 

= -0.41295m^ + 4.9437m + 3.06205 . 

Therefore, it is easy to see that D 2 is nonzero for m > 20. For 2 < m < 19 we again 
refer to Table which completes our proof. 

□ 
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5. Resolving Conjecture 2.10 for small values of m and n 


The main result of this section extends Theorem 14.51 to n < 11, for small m: 


Theorem 5.1 (Resolution of Conjecture 2.10 for small m and n). For m = 2,3,4, 5 
and n = 5, 7,9,11, the network Km,n has the capacity to admit multiple nondegenerate 
positive steady states. 


Proof. We generate our rates and concentrations as in Theorem |3.10| it is straightfor¬ 
ward to check that = 1, A = 1, and e = 0.001 satisfy all necessary hypotheses. Thus, 


we obtain two steady states, r* = (1,1,..., 1) and defined in (131. Then, as in the 
proof of Theorem 3.10 we verify that the determinants det{df{x*)) and det{df{x'^)) 


are nonzero for m = 2,3, 4, 5, which is readily seen from their graphs, which appear in 
Appendix B. (In fact, the graphs strongly suggest that the conjecture holds completely 
for each of these values of n, namely n = 5,7, 9,11, i.e. for m > 5 as well.) □ 


6. Discussion 

As stated in the introduction, deciding whether a chemical reaction network is 
multistationary is not easy in the general case. And even when we can confirm that 
a network is multistationary, there is no general technique to show that it will admit 
multiple nondegenerate steady states. Nonetheless, in this paper we succeeded in this 
task for certain sequestration networks Km,n by using the determininant optimization 
method to obtain closed forms for reaction rates and steady states. 

Our work resolved the n = 3 case of Conjecture |2.10| and we believe that our 
results form an important step toward resolving the full conjecture. Specifically, one 
could use the formulas for rates and steady states given in Theorem |3.10| to analyze 
the general case, or, perhaps easier, the case of some fixed m and general n. Two 
other possible approaches are to (1) find an alternate method to obtain closed forms 
for the steady states of a chemical reaction network, or (2) identify criteria that can 
guarantee that steady states are nondegenerate. 

Expanding on the last idea, our ultimate goal is to develop general techniques 
to assert that steady states of a chemical reaction network are nondegenerate. For 
instance, our analysis of the Jacobian determinants in this work suggest that even 
if the determinant optimization method yields a degenerate steady state, then the 
rate constants can be perturbed slightly so that the degenerate steady state becomes 
nondegenerate (and the other steady state also remains nondegenerate). Is this true 
for any network for which the determinant optimization method applies? If so, then 
this would completely resolve Conjecture |2.10[ and, more generally, this would enable 
us to more readily “lift” multistationarity and thereby enlarge our catalogue of known 
multistationary networks. 
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Appendix A: proofs of Lemmas |4.3| and |4.4[ 


Lemma 6.1 (Lemma 
inequalities for all m > 


4.31. 


The function r 2 {m) 


1.31 
1.31 - 

e —1 


satisfies the following 


m + 1 > r 2 {m) > m . 
Proof. For the upper bound, first observe that 


log 1 -f 


1.31 


m + l 


< log(e^' ) = 1.31 for all m > 0 


m + 1 ^ 

since lima,_^oo(l + -)^ converges to from below for positive values of y. Thus: 


1.31 > log 14- 


1.31 
m -I- 1 


> log 

which implies that 


m + 2.31 
m+l 


m+l 


m+l 


= (m+l) log 


m + 2.31 
m+l 


1.31 
m + l 


> log 


1.31 

e™+i > 


m + 2.31 
m + l 
m + 2.31 
m+l 


=+ (m + 1) — (m + 1) > 1.31. 

and this final inequality implies our desired upper bound: m + l > im^ — = r2(m). 

e m + l 
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Notice that our lower bound is equivalent to the following inequality: 


m + 1.31 131 

- > 


for all m > 2 . 


( 18 ) 


We set a = 1.31, make use of the change of variables z = ^, and then apply log to 
lesired inequalil 

log(l + az) > 


see that our desired inequality (181 is equivalent to: 

az 


a 


+ 1 


1 + z 


for all 2 € (0,1/2) . 


We will show that log(l + az) — > 0. To this end, define &byl — ti = a — 1, and 

notice that l>6>|>(a—1). Next, note that we have the following equalities: 


log(l + az) — 


- ri- — 

Jo \l + t 1 + zj 

= - ^-) + r —^ 

JO (l + j)(l + t) ./bzVl+t 1 + =/ (1 + ^)(1 


+ 1 ) 


(19) 


The second integral in ( |19[ ) is nonnegative (because its integrand is nonnegative), so 
we complete the proof now by showing that the sum of the first and third integrals 


in (|19|) is nonnegative: 
z — t 


f 


(1 + 2)(1 + t) 


dt + 


z — t 
(1 + 2)(1 + t) 


{l — b)z^b —(a —1)^2^ 


(^ + 1)2 


(z + l)^ 


> 0 , 


where the two inequalities come from recalling that 6 < 1, and, respectively, (1 — 6) = 
(a — 1) and & > (a — 1). □ 


Lemma 6.2 (Lemma 


4.41. The function r^{m) = 2 i/Z+sii — 

e ’ ™ + l' -1 


satisfies: 


0.14m > re{m) > 0.13m —0.5 , (20) 

where the upper bound holds for m > 2, and the lower bound holds for m > 20. 


2.1TT1 + 3.41 

Proof. We first prove the upper bound. By the second inequality in (15 I, (e — 


1) > 0 for m > 2. Thus our desired upper bound in (201 is equivalent to the following: 


2.1?n + 3.41 

m(0.14e ™+i - 1.14) > -1.31, 


2.1m + 3.41 

which holds (for positive m) whenever (0.14e ’"+1 — 1.14) > 0. This inequality in 

turn is equivalent to the following (since log is an increasing function): 


2.1m+ 3.41 > (m+l)log(^^ 


(m + l)(2.10) , 


which is true for positive m, so the proof of the upper bound is complete. 

For the lower bound, by clearing the denominator and gathering exponential terms 
on the right-hand side, we see that the desired inequality is equivalent to the following: 

1.13m-1.81 > (0.13m-0.5) . 
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We prove this now. The first inequality below is equivalent to the inequality .00004m > 
—2.536, which is true for positive m: 

1.13m-1.81 > 8.692(0.13m-0.5) 

2 . 1 ( 20 )+ 3.41 2.1771 + 3.41 

> g (2o)+i (0.13m —0.5) > e "*+i (0.13m —0.5) , 

2 . 1777 + 3.41 

and the final inequality holds for m > 20 because e is a decreasing function 

for positive values of m. □ 


Appendix B: graphs for the proof of Thoeorem |5.1| 

Figures [2a}|3a| below present the graphs of the determinant of the Jacobian matrix 
evaluated at the steady states x* and x* (as described in the proof of Theorem 5.11 
for the network Km,n for odd 5 < n < 11 as functions of m. Note that the graphs 
are nonzero for 2 < m < 5, confirming Conjecture |2.10| for odd 5 < n < 11 and those 
values of m. Also, the graphs strongly suggest that the conjecture holds for larger m 
as well. All graphs were made using Desmos Graphing Calculator m- 




(b) det(d/(x#)) 


Figure 2: Km ,5 




Figure 3: Km ,7 
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Figure 4: Km ,9 




Figure 5: Km ,11 
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